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ABSTRACT 

Dynamics  of  the  circulation  in  the  Sea  of  Marmara  are 
investigated  with  a  time  dependent,  three  dimensional 
numerical  model.   The  empirically  inferred  hydrologic  regimes 
of  the  sea  and  connective  straits  are  discussed. 

A  baroclinic  model  based  on  the  primitive  equations  is 
solved  by  direct  integration  of  an  initial  value  problem. 
The  circulation  in  the  Sea  is  driven  by  surface  forces  that 
simulate  wind  stress  and  horizontal  pressure  gradient  forces 
related  to  internal  stratification. 

The  predicted  density  field,  in  sigma— t  units, is  compared 
with  data.   Detailed  three  dimensional  horizontal  velocity 
patterns  and  vertical  velocity  patterns  in  horizontal  planes 
are  given. 

Bottom  friction,  irregular  bottom  topography,  non- linear 
terms  in  the  momentum  equation  and  vertical  mean  part  of  the 
horizontal  velocity  have  been  omitted.   For  simplicity 
density  is  predicted  in  place  of  temperature  and  salinity. 
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I.   INTRODUCTION 

There  is  an  intense  growing  interest  in  recent  years  in 
the  development  of  simulation  of  hydrological  conditions  in 
coastal  waters  and  smaller  adjacent  seas.   Advancement  in 
modern  technology  requires  much  greater  accuracy  than  can  be 
obtained  by  the  use  of  classical  instrumentation  and  measure- 
ment techniques  only.   Due  to  the  discontinuous  nature  of 
such  observations  these  measurements  can  provide  the  basis 
for  analysis  of  the  current  system,  temperature,  salinity 
and  other  variables  on  a  climatological  basis.   But  climato— 
logical  conditions  are  not  accurate  enough  in  coastal  waters 
for  uses  such  as  construction,  navigation,  tracing  of  pollu- 
tants, mining  effluent,  oil  spills,  search  and  rescue  opera- 
tions, agriculture,  etc.   It  is  also  true  that  sampling  the 
ocean  adequately  whether  for  exploratory  or  research  purpose 
is  a  major  problem.   Great  short— term  variability  of  condi- 
tions in  coastal  waters  makes  it  even  more  difficult  to  make 
adequate  representative  measurements. 

As  a  result,  in  addition  to  the  classical  methods  of 
analyses  of  current  systems,  determination  of  salinity  and 
temperature  distributions  etc.,  new  simulation  techniques 
which  were  applied  long  before  synoptic  numerical  weather 
prediction  received  much  more  attention  (T.  Levastu,  S. 
Larson  and  K.  Rabe,  1976) . 
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The  objective  of  this  study  is  to  develop  a  regional 
circulation  model  of  the  Sea  of  Marmara.   The  circulation 
and  density  distribution,  in  sigma— t  units,  in  that  sea  are 
simulated  by  a  nine- level  numerical  model  with  an  hori- 
zontal flat  bottom  at  1000  m  depth.   The  motion  is  driven 
by  prescribed  wind  stress  and  horizontal  pressure  gradient 
forces  due  to  internal  stratification. 
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II.   PHYSICAL  GEOGRAPHY  AND  HYDROLOGY  OF 
THE  SEA  OF  MARMARA 

The  Sea  of  Marmara  is  a  small  sea  located  between  Asia 
Minor  and  Europe  and  occupies  the  main  part  of  the  Turkish 
straits  system  (210  km) .   The  Strait  of  Canakkale   (Darda- 
nelles) (60  km)  connects  the  Sea  of  Marmara  to  the  Mediter- 
ranean and  through  the  Strait  of  Istanbul  (Bosporus)  (30  km) 
it  is  connected  to  the  Black  Sea  at  the  North.  (Figure  1) 

The  structure  of  the  basin  which  is  small  in  area  (about 

2 
11000  km  )  is  associated  with  the  Anatolicin  fault  which 

lies  along  its  bottom  to  the  North.   The  trough  passing  along 
the  Northern  deep  slope  consists  of  three  depressions  (max. 
1380  m) .   These  three  basins  are  separated  by  low  connecting 
sills.   A  shallow  trough  lies  at  the  foot  of  the  continental 
slope  which  is  well  pronounced  only  in  the  easternmost  part 
of  the  sea.   At  the  South  it  becomes  shallower  having  depths 
of  the  order  60—80  m.   Minimum  depths  at  the  Strait  of 
Canakkale   and  Strait  of  Istanbul  are  57  and  37  m;  they  have 
mean  widths  4.5  km  and  0.7  km  and  average  lengths  60  and  30 
km,  respectively. 

The  climate  of  the  Sea  of  Marmara  is  influenced  by  the 
regional  atmospheric  circulation  system  of  the  Eastern  Medi- 
terranean, the  Black  Sea  and  Asia  Minor  Peninsula.   During 
the  summer,  the  Eastern  Mediterranean  basin  is  under  the 
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influence  of  a  stable  wind  system  which  is  called  "Etesian" 
wind.   Etesians  are  monsoon— like  and  associated  with  a  quasi- 
stationary  cyclonic  circulation  above  the  Asia  Minor  Penin- 
sula during  summer.   The  active  period  of  the  Etesian  winds 
falls  almost  entirely  within  the  summer  season  (approximate- 
ly from  May  to  October) .   There  are  two  maximums,  which  are 
pronounced  during  July  and  August.   These  Etesian  winds, 
which  are  dry  winds  and  generally  northerly  or  northeasterly 
from  the  Black  Sea,  have  diurnal  variability. 

Southerly  winds  bring  warm  and  humid  air  to  the  area 
during  winter.   When  these  southerly  winds  are  sufficiently 
strong,  they  produce  storm  waves  and  surges  along  the  coast 
of  Marmara  and  at  the  entrance  to  the  Strait  of  Istanbul 
(Gunnerson  and  Ozturgut,  1974) . 

The  hydrological  regime  of  the  Sea  of  Marmara  is  largely 
determined  by  the  water  exchange  between  the  Mediterranean 
and  Black  Seas  through  the  connective  straits,  Strait  of 
Canakkale,  and  Strait  of  Istanbul.   Problems  of  water  ex- 
change, vertical  mixing  and  hydrological  characteristics  of 
the  currents  in  these  straits  and  in  the  Sea  of  Marmara  have 
been  the  subject  of  many  extensive  studies  conducted  by  the 
Turkish  Navy  Hydrographic  Office  and  by  the  University  of 
Istanbul.   In  addition  to  these  fundamental  and  frequent 
studies  there  were  some  investigations  conducted  by  German, 
Russian  and  English  investigators  at  the  end  of  the  19th 
and  beginning  of  the  20th  century.   Recent  studies  and  obser- 
vations, if  connected  to  old  data  and  observations,  can  give 
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a  better  understanding  of  the  complete  regional  oceanographic 
condition. 

These  studies  are  mainly  concentrated  on  the  Strait  of 
Istanbul  and  the  vicinity  of  this  strait.   This  is  due  to  the 
importance  of  this  strait  as  a  link  in  the  connection  between 
the  Black  Sea  and  the  Mediterranean  Sea  through  the  Sea  of 
Marmara.   This  makes  it  an  important  factor  controlling  Black 
Sea  chemistry,  biology,  sedimentation,  etc. 

1  V 

-*>  Based  on  the  first  detailed  survey  of  Makarov  and  Merz 

1917-18,  in  Moller  (1928),  and  subsequent  studies  of  Ullyott 

-*   and  Ilgaz  (1946),  Pektas  (1956),  A.  K.  Bogdanova  (1961,  1965), 
C.  G.  Gunnerson  and  E.  Ozturgut  (1974)  and  others,  the  main 

-y      properties  of  this  straits'  system  and  its  vicinities  are 
fairly  well  understood. 

One  of  the  main  characteristics  of  these  straits  is  a 
two— layer  structure  of  the  current,  which  is  determined  by 
the  differences  between  the  relatively  fresh  water  of  the 
Black  Sea  and  saline  water  from  the  Mediterranean  which  occu- 
pies the  depths  of  the  Sea  of  Marmara  and  forms  the  bottom 
water  of  two  connective  straits  (Figure  2).   Other  external 
factors  such  as  wind,  water  level  drop  along  the  straits, 
water  balance,  topographic  influences  and  the  depth  and  width 
of  the  straits  also  must  be  considered.   The  scale  of  the 
motion  and  complexity  of  these  factors  make  it  difficult  to 
determine  the  current  system  mathematically. 

Another  main  characteristic  of  these  straits  is  the  two- 
layer  stratification  of  water  which  is  associated  with  the 
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two— layer  current  system.   This  stratification  is  also  a 
consequence  of  the  water  mass  differences  at  the  two  ends 
of  the  straits. 

Detailed  observations  of  temperature  and  salinity  dis- 
tributions indicate  certain  features  of  the  diffusion  of 
the  Mediterranean  water  into  the  Black  Sea.   This  water 
flows  North  through  the  Strait  of  Istanbul  from  the  Sea  of 
Marmara.   It  is  shown  by  Bogdanova  (1961,  1965)  that  this 
underflow  of  Mediterranean  water  into  the  Black  Sea  is 
present  throughout  the  year  and  descends  to  significant 
depths  in  the  Black  Sea. 

This  water  exchange  through  the  straits  system  fluctu- 
ates according  to  seasons,  years  and  periods.   This  is  very 
important  for  the  hydrologic  regimes  of  the  Sea  of  Marmara 
and  the  Black  Sea.   Water  level  differences,  given  with  the 
following  expression,  are  mainly  due  to  excess  precipita- 
tion over  evaporation  and  river  runoff  in  the  Black  Sea. 


Ah  =  ZB  -  ZM  (1.1) 


Ah  :  sea  level  differential  between 
Black  Sea  and  Mediterranean 


ZR  :  water  level  of  Black  Sea 

Z.,  :  water  level  of  Mediterranean 
M 


Sea  level  differential  (Ah)  is  always  positive  and  fluctu- 
ates with  time.   It  can  be  expressed  in  two  components 

Ah  =  Ah  +  (Ah)  »  (1.2) 
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where 


Ah  :  average  sea  level  differential 

(Ah) ' :  fluctuating  part  of  sea  level 
differential  [(Ah)'  «  f(t)] 


According  to  Bogdanova  (1965)  Ah  has  a  value  42  cm 
based  on  the  long  term  observations.   The  fluctuating  part 
of  the  sea  level  differential  and  total  level  differential 
are  given  in  Table  I. 

TABLE  I. 

Monthly  average  fluctuation  sea  level  differential 
and  sea  level  differential. 1 


Month 

I 

II 

III 

IV 

V 

VI 

VII 

VIII 

IX 

X 

XI 

XII 

(Ah)  '  (cm) 

1 

0 

3 

7 

11 

15 

6 

0 

-5 

-7 

-7 

-4 

Ah  (cm) 

43 

42 

45 

49 

53 

57 

48 

42 

37 

35 

35 

38 

(  Based  on  the  data  given  by  Bogdanova,  1965) 


The  maximum  change  is 


-^  (Ah)   v  =  ^r  (Ah')   v  =  -  9  cm/month, 
3t     max    9t      max 


observed  during  approximately  June.   This  falls  approximately 
in  a  period  when  continental  runoff  begins  to  decrease.   The 
fluctuating  sea  level  differential  increases  from  March  to 
July.   As  a  consequence  during  this  period  the  upper  current 
becomes  stronger  and  the  lower  current  becomes  weaker.   From 
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August  to  November  sea  level  fluctuation  decreases  and  be- 
gins to  increase  again  from  August  to  February.   The  upper 
current  of  the  Strait  of  Istanbul  becomes  weaker  and  the 
lower  current  becomes  stronger. 

According  to  Bogdanova's  (1965)  observations  the  maxi- 
mum Ah  value  is  observed  during  June  and  has  a  value  57 
cm;  the  minimum  Ah  value  35  cm  is  observed  during  October 
and  November. 

The  average  sea  level  differential  between  the  Black 
Sea  and  the  Sea  of  Marmara  is  about  3  5  cm  which  is  six 
times  Moller's  (192  8)  estimate  (Gunnerson  and  Ozturgut, 
1974) . 

There  is  a  characteristic  stratification  with  a  sharp 
density  transition  layer  over  the  whole  area.   The  sharpness 
is  more  pronounced  at  the  North  at  the  Strait  of  Istanbul. 
A  cross  section  of  sigma— t  values  for  the  Strait  of  Istanbul 
is  shown  in  Figure  3.   Fresh  warm  water  of  the  Black  Sea 
lies  in  a  surface  layer  15—20  m  thick  at  Istanbul.   This 
water  which  is  mainly  river  runoff  and  excess  precipitation 
overrides  the  Mediterranean  water  which  is  cold  and  dense. 
The  wedge  form  of  the  upper  water  shows  clearly  in  this 
strait,  but  this  is  not  true  for  the  Strait  of  Canakkale. 
On  the  other  hand,  the  wedge  form  of  the  lower  water  is  com- 
paratively stronger  and  more  pronounced  in  the  Strait  of 
Canakkale  due  to  the  downward  sloping  of  the  sea  bed  toward 
the  North  in  the  Strait  of  Istanbul.   Average  sigma— t  values 
for  surface  and  bottom  water  are  13.5  and  27.5.  (Defant,  1961) 


19 


p 

XI 

c 

(C 

+J 

CO 

H 

«4-l 

0 

-M 

•H 

id 

n 

+j 

w 

<w 

0 

C 

• 

0 

*-^ 

-H 

CO 

-P 

4-1 

U 

c 

<D 

•H 

CO 

0 
o 

CO 

i-ij 

CO 

Q 

0 

1 

u 

U 

o 

c 

-p 

CD 

1 

0) 

<d 

5 

£  +> 

Cn  O 

•H  XJ 

w 

"— * 

m 

•H 

fa 


-i 1 r 

°(i/\i)  h    i  °d    3    a  * 


o 


20 


According  to  Gunnerson  and  Ozturgut  (1974) ,  the  average 
salinity  of  Black  Sea  waters  at  a  station  9  km  from  the 
Northern  entrance  of  the  Strait  of  Istanbul  was  approximate- 
ly 17.5  °/oo  for  an  observation  period  July  1966  to  December 
1967.   Lower  values  were  observed  (16—17  /oo)  from  July  17 
to  August  23,  1967.   This  lower  surface  salinity  reflects 
maximum  discharge  of  rivers  into  the  Black  Sea.   Higher 
average  surface  salinities  (17.5  —  19.0   /oo)  were  observed 
for  the  rest  of  the  year  and  extreme  values  (as  high  as 
25  /oo)  during  winter  months.   These  were  due  to  the  favor- 
able conditions  for  Mediterranean  water  to  flow  northward. 
These  conditions  include  small  differences  between  the  levels 
of  the  Black  Sea  and  the  Mediterranean  and  southerly  or 
southwesterly  winds  over  the  Sea  of  Marmara  and  Strait  of 
Istanbul.   At  the  midpoint  the  average  salinity  in  the  lower 
layer  is  approximately  38.5   /oo  and  in  the  surface  layer 
17.5   /oo.   A  salinity  cross  section  is  given  in  Figure  4. 

Warm  and  fresh  Black  Sea  water  can  be  identified  easily 
from  the  temperature  cross— section  which  is  given  in  Figure 
5.   Below  this  warm  fresh  water  of  the  Black  Sea  cold  inter- 
mediate water  of  Sea  of  Marmara  can  be  observed.   Mediter- 
ranean water  occupies  the  deep  bottom  layer  of  the  Strait  of 
Istanbul. 

The  dynamics  of  the  processes  occurring  in  these  two 
straits  is  determined  by  the  wind  stress  over  the  straits 
and  the  drop  of  the  sea  levels  and  densities  at  the  end  of 
the  straits.   The  main  driving  forces  are  gradient  of  water 
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level,  the rmo ha line  forcing,  wind  force  exerted  on  the  sur- 
face and  coriolis  force. 

According  to  observations  the  boundary  between  the  water 
masses  generally  does  not  coincide  with  the  boundary  between 
currents.   The  boundary  between  water  masses  has  a  greater 
slope  than  the  boundary  between  currents  (Defant,  1961) . 
This  is  probably  due  to  diffusion  processes  occurring  at  the 
current  boundary  and  to  the  nature  of  the  current  (Figure  2) . 

The  average  Black  Sea  inflow  and  outflow  through  the 

Straits  of  Istanbul  according  to  the  analysis  of  Moller  (1928) 

3 
are  6,100  and  12,600  m  /sec,  respectively  (Gunnerson  and 

Ozturgut,  1974) .   Velocity  is  greatest  at  the  sea  surface 

and  decreases  rapidly  with  depth  and  laterally;  it  increases 

from  North  to  South.   Under  normal  conditions  it  is  40—50 

cm/sec  at  the  entrance  of  the  strait  and  150  cm/sec  at  the 

South  end.   The  lower  current  is  strongest  in  the  central 

part  of  the  lower  water  which  is  at  about  16  m  from  the 

bottom  in  the  Strait  of  Istanbul  and  45  m  in  the  Strait  of 

Canakkale.    These  velocities  are  100—150  cm/sec  and  25  to 

10  cm/sec  respectively  (Defant,  1961) . 

In  general  in  the  Sea  of  Marmara  there  is  Black  Sea 

water  at  the  surface  and  Mediterranean  water  at  the  bottom. 

This  is  the  normal  state.   Except  for  times  when  very  strong 

southerly  winds  bring  bottom  water  to  the  surface  and  cause 

wind  mixing  these  two  main  water  types  are  both  present  with 

a  very  thin  intermediate  water  layer.   The  Black  Sea  water 

extends  to  a  depth  of  15—20  m  in  the  vicinity  of  the  Strait 
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of  Istanbul  and  becomes  shallower  farther  south.   There  is 
intermediate  water  between  these  two  water  types  which  ex- 
tends at  most  to  about  50—60  m.   Below  this  the  basin  is 
filled  with  Mediterranean  water  as  shown  in  Figures  6,7  and 
8.   This  water  has  a  temperature  14.9°C,  salinity  38.55   /oo 
and  sigma— t  28.75,  approximately  in  June. 

Surface  water  has  a  large  annual  variability.   On  the 
other  hand  bottom  water  does  not.   Due  to  the  described 
oceanographic  conditions  bottom  water  is  almost  isolated  from 
surface  processes  and  has  a  large  time— scale  of  variability 
compared  to  surface  water.   Average  inflow  and  outflow  rates 
at  the  straits  give  approximately  30  years  for  "turnover" 
time  or  "flushing"  time  for  bottom  water  and  200  days  for 
surface  water. 

Based  on  the  same  rates  a  typical  value  for  the  vertical 
mean  current  |u|  is  10    cm  sec   .   A  characteristic  magni- 
tude for  horizontal  velocity  |u|  using  the  thermal  wind 
relation  is  4  cm  sec   .   Due  to  the  large  isolated  reservoir 
of  Mediterranean  water  of  almost  uniform  characteristics  in 
the  lower  part  of  the  basin,  the  thermohaline  circulation 
is  confined  to  the  upper  50—60  m. 
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III.   SELECTION  OF  THE  NUMERICAL  MODEL 

Regional  circulation  models  may  be  separated  into  two 
categories.   One  type  of  regional  circulation  model  is  a 
diagnostic  density  model.   This  approach  is  based  on  an  ob- 
served density  field  and  wind  stress  distribution.   These 
quantities,  however,  are  not  observed  continuously  and  in 
an  evenly  distributed  way.   But  based  on  average  values  of 
these  quantities  a  velocity  field  is  calculated  diagnosti— 
cally.   A  diagnostic  model,  which  seems  at  first  easier, 
has  drawbacks  mainly  of  two  types.   First,  it  requires  tak- 
ing into  account  actual  configuration  of  the  basin  geometry 
and,  second,  it  requires  sufficient  representative  observa- 
tions.  These  two  factors  make  this  approach  complicated 
and  expensive. 

A  second  type  of  approach  to  study  regional  circulation 
models  is  to  use  a  predictive  model  with  an  idealized  topo- 
graphy and  basin  area.   But  this  becomes  more  complex  by  the 
inclusion  of  more  dynamical  processes. 

At  first  glance  it  seems  that  diagnostic  models  are  a 
better  way  of  simulation  of  dynamical  processes.   But  these 
observations  on  which  such  models  depend  are  generally 
gathered  by  land  based  stations  (e.g.,  wind  data)  and  extra- 
polated to  the  area  of  concern.   Furthermore,  using  average 
values  brings  up  the  question  of  how  representative  the 
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result  is.   Further,  it  is  difficult  to  get  a  continuous 
view  of  the  dynamical  process  and  long— term  quantative 
analyses. 

On  the  other  hand,  although  prognostic  models  are  complex 
due  to  the  dynamical  processes  involved,  they  are  less  expen- 
sive and  can  give  a  continuous  view  of  variables,  which  is 
generally  based  on  an  initial  density  field.   At  least  as  a 
first  approximation  it  seems  reasonable  and  economical  to 
study  a  regional  circulation  simulation.   Present  data  are 
still  used  to  estimate  the  initial  density  field  and  wind 
stress  distribution  in  this  kind  of  model. 

Due  to  the  uneven  distribution  of  density  and  wind  stress 
observations  it  was  mandatory  to  use  the  prognostic  approach 
at  the  beginning.   The  procedure  adopted  is  to  use  present 
data  as  efficiently  as  possible  to  estimate  a  realistic  ini- 
tial density  field  and  a  constant  wind  stress  distribution 
throughout  the  integration  period. 
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IV.   DESCRIPTION  OF  THE  DYNAMIC  MODEL 

A.   BASIC  EQUATIONS  OF  THE  MODEL 

Basic  assumptions  which  lead  to  important  simplifications 
are:  (1)  the  hydrostatic  assumption,  (2)  the  Boussinesq  assump- 
tion with  respect  to  density  variations  (density  variation  is 
neglected  except  where  it  appears  as  a  coefficient  of  the 
gravitational  constant) ,  and  (3)  an  implicit  treatment  of  the 
mixing  process  by  eddy  diffusion.   In  addition,  the  non— linear 
terms  are  neglected  in  the  momentum  equation. 

With  these,  the  sea  is  assumed  incompressible,  and  density 
is  replaced  by  a  constant  value  everywhere  except  in  the  hydro- 
static equation. 

The  governing  equations  based  on  these  assumptions  are: 


|u  =  _  8  (i,  +  a  v2u  +  K  ^  +  fv 

8t      3x  p      M       v  ~  2 
Ko  8z 


(4.1) 


9v               3   /  p  \         ah2           „  3    v         ,.                                               ,,    ,, 

Tt  =  "  3y(p-°    +AMVv  +   Kv~2-fu  (4-2) 

J         O  dZ 

H--(p-p0)g  (4.3) 

iH  +  iZ  +  i^=0  (4    4) 

H-  V2p  +  KH+  *olp)  (4-5) 

dZ 
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where  the  symbols  have  the  following  meanings: 

t  time 

x  horizontal  coordinate  measures  positive  eastward 

y  horizontal  coordinate  measures  positive  northward 

z  vertical  coordinate  measures  positive  upward  with 
an  origin  at  mean  water  level 

u  eastward  velocity  component 

v  northward  velocity  component 

w  vertical  velocity  component 

p  density 

p  reference  constant  density 

P  departure,  from  reference  pressure  which  is 

P  pressure  at  reference  density  (=  —  Pogz) 

g  acceleration  of  gravity 

f  coriolis  parameter  (=  2ft  sin<j>) 

ft  angular  speed  of  earth  rotation 

<j>  latitude 


AM 


coefficient  of  horizontal  diffusion  for  momentum 
(constant) 


K       coefficient  of  vertical  diffusion  for  momentum 
(constant) 

A_^      coefficient  of  horizontal  diffusion  for  density 
(constant) 

K       coefficient  of  vertical  diffusion  for  density 
(constant) 

iP  material  derivative  [=  IP  +   u  ii}+  v  ^  +   w  |i}  ] 

dt  3t      8x      8y      dz 

6       introduces  the  effect  of  convective  process  when 
the  stratification  is  unstable 
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The  sea  is  driven  by  wind  and  thermohaline  forcing,  these 
effects  being  introduced  with  surface  boundary  conditions  in 
the  momentum  and  density  equations,  respectively.   The  bound- 
ary conditions  specified  at  the  sea  surface  are 


KH  -o 


K  lH=o 

Kv  8z    U 


y       /     z  =  0  (4.6) 

T. 


v  az  "  po 

w  =  0 


Heating  or  cooling  at  the  sea  surface  is  neglected  over 
the  period  of  integration.   A  linear  north— south  density  gra- 
dient is  specified  in  the  uppermost  level  in  the  model  and 
held  constant  during  calculations.   The  x  component  of  sur- 
face stress  is  taken  to  be  zero  and  calculations  are  carried 
out  with  a  uniform  meridional  wind  stress. 

The  boundary  conditions  specified  at  the  flat  horizontal 
bottom  are 

K  l£  =  0 

(    z  =  -H  (4.7) 

u=v=w=  0  1 
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The  boundary  conditions  at  the  side  walls  of  the  basin 
are 


Ah  |*  =  u  =  v  =  0       J     y  =  0,  Ly 


3y 


AR  |g  =  u  =  v  =  0 


I 


where      Ly  :  width  of  the  basin  (  50  km) 
Lx  :  length  of  the  basin (200  km) 

Boundary  conditions  at  the  straits  are  based  on  specified 
density  and  velocity  fields  which  were  taken  from  available 
data  and  annual  average  inflow  and  outflow  rates  respectively. 

B.   DESCRIPTION  OF  THE  SOLUTION  TECHNIQUES 

Since  the  model  sea  has  a  horizontal  flat  bottom  and 
"rigid  lid"  approximation,  vertical  velocity  vanishes  at  the 
surface  and  bottom.   Following  the  method  used  by  Bryan  and 
Cox  (1967-1968a-b) ,  Haney  (1974)  and  others,  the  ocean 
surface  was  considered  a  balanced  surface  at  which  w  =  0. 
In  this  approach  the  height  of  the  sea  surface  is  not  obtained 
from  the  continuity  equation;  therefore,  external  gravity 
waves  are  eliminated.   This  also  removes  the  vertical  mean 
divergence  from  all  scales  of  motion.   Elimination  of  external 
gravity  waves  allows  the  use  of  a  longer  time  step.   In  the 
case  of  long— term  integration  this  is  an  important  gain  com- 
pared to  the  loss  of  the  prediction  for  height  of  the  sea 
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surface.   Filtering  external  gravity  waves  has  a  negligible 
effect  on  the  density— driven  part  of  the  current. 

Due  to  the  conditions  on  vertical  velocity,  integration 
of  the  continuity  equation  from  bottom  to  surface  requires 
the  following  condition: 


lH  +  IX  =  o  (4  9) 

8x  +  By   u  l4,y; 


where  the  meaning  of  the  overbar  is 


1  r° 
(  )  vertical  mean  [=  —   /   (  )  dz] 


and 


-H 


(  ) ' indicates  departure  from  vertical  mean 
of  any  quantity. 


Also  for  some  basic  integral  constraints  this  is  a  required 
relation  to  be  consistent  with  boundary  conditions. 

The  "rigid  lid"  approximation,  and  (4.9) ,  are  accomplished 
by  vertically  integrating  equations  (4.1),  (4.2)  and  subtract- 
ing the  result  from  equations  (4.1)  and  (4.2)  respectively. 
Using  present  boundary  conditions,  the  new  sets  of  equations 
are 

9u»  8P'         _    n2    ,         _      82u'         .    , 

TT-r—   = r—   +    A..V    u"     +    K      «•      +    fv' 

dt  p    8x  M  v    -    2 

O  dZ 


(4.10) 


y 

vr       = 5-7    +    A..V     V'     +     K        x- rr    —    fu1 

at  p    8y         Tl  v   ,    2  pH 


o   J  "  v    3z"  Mo* 
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where 

u1  eastward  shear  velocity  component 

v'  northward  shear  velocity  component 

P"  departure  of  P  from  its  vertical  average 

[=  J     (p-po)gdn  -  |  J     (  f     (p-po)gdn)dz]       (4.11) 
z  —Hz 

The  relation  between  shear  current  and  total  current  for 
horizontal  components  is 


u* (x,y,z,t)  =  u(x,y,z,t)  -  u(x,y,t) 


v' (x,y,z,t)  =  v(x,y,z,t)  -v(x,y,t) 


(4.12) 


since 

|u|  >  >  |u|   and   |v|  >  >  |v|  (4.13) 

equation  (4.12)  takes  the  form 

u' (x,y,z,t)  =  u(x,y,z,t) 


(4.14) 


v' (x,y,z,t)  =  v(x,y,z,t) 


Thus  total  u,  v  component  of  horizontal  velocity  can  be 
represented  by  vertical  shear  currents  u',  v'.   In  the 
following  part  primes  of  u'  and  v'  are  dropped. 


36 


The  equations  are  non— dimensionalized  by  making  the  fol- 
lowing substitutions  for  new  non-dimensional  variables 

(x,y)  =  (x,y)/L 
^  (z)    =  (z)/h 
(u,v)  =  (u,v)/V, 


where 


(w) 

(vO/V^ 

~L~ 

(P,P')= 

(P,PM 
PofLVl 

(t) 

(t) 

L 
Vl 

(T)     = 

(T) 

TM 

(a) 

(a) 

(aN"aS) 

(4.15) 


V,     scale  velocity  associated  with  density  driven 
current 

h     scale  depth,  taken  equal  to  a  characteristic 
thermocline  depth 

L     length  scale 

xM    wind  stress  [=  MAX(x  (x,y) ,  xy(x,y)] 

3 
a  sigma— t  [=  (p— 1)  x  10  ] 

a  surface  sigma— t  at  northern  part  of  the  basin 

aq    surface  sigma— t  at  southern  part  of  the  basin 
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A  scale  velocity  associated  with  the  density  driven  cur- 
rent is  defined  by  using  geostrophic  and  hydrostatic  scaling 
which  were  introduced  by  Bryan  and  Cox  (1968a) .   When  the  geo- 
strophic relation  is  differentiated  with  respect  to  z  and  the 
hydrostatic  equation  is  used,  the  result  is 


o 


A  scale  velocity  connected   with  the  density  distribution 
may  be  defined  by 


V,  =2_A£h  (4.17, 


where 


Ap    north— south  surface  density  differences 
(=  PN  -  PS) 


A  scale  velocity  for  wind  driven  current  is  defined  by 


V0   =  -   "  (4.18) 

2     P0fhLr 


where 


x..    wind  stress 
M 

L     relative  characteristic  length  scale 
for  wind  stress  [=  L/R]  and  L  is  a 
length  scale  for  wind  stress  and 
R  radius  of  the  earth 
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Substitution  of  these  non-dimensional  variables  into 
(4.10),  (4.11),  (4.4)  and  (4.5)  with  some  rearrangement, 
gives  these  equations  in  terms  of  non-dimensional  variables 
of  the  form 

3u      1  3P*    EH  n2    EV  92u    1  ,,  1Q, 

7t  "  ~  Ito  ?x   +Ro-Vu+Ro-7-2  +  Ro-V  (4'19) 

dZ 

Sz  =  -  jl  3P'  +  fs  y2v  +  fx g2v  -  ^  _i_  t  y 

dz  Ro  8y    Ro       Ro  ~7~2        Ro  L  H   s 

dz        r  r 

-Au  (4.20) 

o  o    o 

P«  =  J    (a-aQ) dn  ~  k   f  {   f   (a-a0>dn)dz         (4.21) 
z  -Hz, 

w  -  J       V-  V  dn  (4.22) 


8a     3a     9a     3a    1  „2     ,  K  .  EV  82a 

Tt  +  u^  +  v^  +  W^=p¥Va+  (!T)  ro"  ~ 2 

V       dZ 

+  6(a)  (4.23) 


where 


Ro   Rossby  number  which  shows  the  relative  importance  of 
local  time  change  term  in  the  equation  of  motion  with 
respect  to  that  of  the  coriolis  term       V, 

I"  IT1 

E   horizontal  Ekman  number  which  shows  relative  impor- 
tance of  lateral  diffusion  and  coriolis  terms 

AM 

fir 

E-.  vertical  Ekman  number,  gives  the  ratio  of  vertical 
diffusion  of  momentum  to  the  coriolis  term 

fh 
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Pe  Peclet  number,  shows  importance  of  advection  compared 
to  diffusion  of  density 

V  L 

V   relative  velocity,  ratio  of  barotropic  velocity  to 
baroclinic  velocity 

v2 

H   relative  depth,  ratio  of  basin  depth  to  thermocline 
r  depth 

L   hJ 


The  model  is  governed  by  these  six  parameters  assuming  a 
value  for  L   is  defined.   All  the  parameters  of  the  resultant 
run  and  constants  of  the  model  are  given  in  table  2. 

Equations  (4.19)  —  (4.23)  are  solved  for  variables  u,  v, 
P',  w  and  a.   The  variables  u,  v,  and  a  are  predicted  from 
(4.19),  (4.20)  and  (4.23).   P'  and  w  are  calculated  diagnosti- 
cally  from  equations  (4.21)  and  (4.22),  respectively. 
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TABLE  2.   PARAMETERS  AND  CONSTANTS 


—3 

Ro      4.255  x  10  Rossby  number 

—2 

E„      0.532  x  10  Horizontal  Ekman  number 

H  —2 

E       0.851  x  10  Vertical  Ekman  number 

Pe      10.0  Peclet  number    -?— 

H       50.0  Relative  depth 

V       0.1  Relative  velocity 

r  5  . 

Ax      10  x  10   cm  Zonal  grid  spacing 

5  .... 

Ay      5  x  10   cm  Meridional  grid  spacing 


r 


10.0 

50.0 

0.1 

10  x 

105  cm 

5  x  105  cm 

8  minute 

28.9 

day 

40° 

2  it  day 

6.37 

x  10   cm 

At      8  minute  Time  step 

t       28.9  day  Time  scale 

<jj       40°  Reference  latitude 

Q  2tt  day  Rotation  rate  of  the  earth 

g 

R       6.37  x  10   cm  Radius  of  the  earth 

—2 

g       1000  cm  sec  Gravity  acceleration 

—3 

p       1.012  gm  cm  Reference  density 

H       1  x  10   cm  Depth  of  the  sea 

L       10  Curl  factor  or  relative  length 


scale  for  wind  stress 


2    —1 
K       1.6  cm   sec  Eddy  diffusivity 

2    —1 
K       3 . 2  cm  sec  Eddy  viscosity 

v  —1 

V-,       4  cm  sec  Characteristic  baroclinic 

velocity 

Characteristic  length  scale 

AM      0.5  x  10°  cm'  sec  """   Lateral  eddy  viscosity 

7    2—1 
A„      0.4  x  10   cm   sec    Lateral  eddy  diffusion  coefficient 


4  cm  sec 

7 
1  x  10   cm 

0.5  x  108  cm2 

-1 
sec 

0.4  x  107  cm2 

-1 
sec 
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V.   DESCRIPTION  OF  THE  NUMERICAL  MODEL 

A.   SPACE  AND  TIME  DIFFERENCING  TECHNIQUES 

The  space  and  time  differencing  schemes  are  similar  to 
the  schemes  used  by  Haney  (1974)   for  a  general  circulation 
model. 

In  the  integration  of  the  primitive  equations  a  staggered 
grid  and  centered  space  differencing  scheme  are  utilized. 
The  main  concern  in  choosing  the  staggered  grid  was  avoidance 
of  a  computational  mode  which  appears  when  a  centered  space 
differencing  scheme  is  used  with  an  unstaggered  grid.   This 
computational  mode  in  space  is  not  present  for  a  staggered 
grid.   The  staggered  grid  also  has  advantage  in  saving  comput- 
ing time  due  to  the  presence  of  variables  at  alternate  grid 
points. 

A  leapfrog  scheme  is  used  for  time  integration.  The 
starting  scheme  is  a  combination  of  forward  and  leapfrog 
scheme  as  shown  below 


Leapfrog  ^At 


n=0     _  ,  ,  Afc  2  n=l 

Forward  %At 


(•=:  time   step   forward)    +    (•=■  time   step    leapfrog)    -*•  Time   step   1 
where 


n      represents   a   time   step 
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In  most  circulation  models  Matsuno  scheme  (Matsuno,  1966) 
is  used  as  a  starting  scheme  and  to  remove  the  solution 
separation  which  is  present  when  a  leapfrog  scheme  utilized 
exclusively  for  time  integration.   Since  time  integration  in 
this  model  is  not  long  enough  to  cause  solution  separation 
a  cheme  such  as  Matsuno  is  not  used  here. 

Each  term  in  the  momentum  and  sigma— t  equations  is  evalua- 
ted at  a  different  time  step  to  ensure  linear  computational 
stability. 

Neglecting  coefficients,  the  time  differencing  of  the 
equation  of  motion  corresponding  to  (4.20)  is  given  by 


vn+l)  =  v(n-l)  +  2At[pT(n)  +  FT(n-D } 


-   2At  f  [CI  u(n+1)  +  C2  u(n  1}]  (5.1) 


where 

PT    pressure  term 

FT    friction  term,  which  is  composed  of  lateral 

diffusion  of  momentum  and  vertical  eddy  stress 

A  trapezoidal  implicit  scheme  is  utilized  for  the  coriolis 

term.   Then,  to  satisfy  consistency  and  linear  computational 

stability,  coefficients  CI  and  C2  must  satisfy  the  relation. 

CI  +  C2  =  1 

j  £  cl  1  1  <5-2> 

The  values  Cl  =  0.55  and  C2  =  0.45  are  commonly  accepted 
for  these  constants  which  produce  very  slight  damping  of 
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inertial  oscillations.   Equation  (5.1)  and  an  analogous  equa- 
tion for  u      are  solved  simultaneously  for  the  two  unknowns 
u     ,  v     .   The  time  differencing  form  of  the  sigma— t 
equation  has  the  form 


a 


(n+1)    =  a(n-l)    +   2At[ADV(n)+  HD<n-l>    +  VD(n-l)] 


+    8c(on+1)  (5.3) 


where 

ADV  advection  term 

HD  horizontal  diffusion  term 

VD  vertical  diffusion  term 

The  leapfrog  scheme  is  conditionally  stable  and  the  time 
step  must  satisfy  the  relation  (for  two  dimensional  wave 
propagation) 

C  —  <  —  (5  4) 

CAx±^  (5-4) 


where 


maximum  phase  speed  of  the  waves  present 
in  the  system. 


External  gravity  waves  are  removed  and  inertial  oscilla- 
tions are  rendered  neutral  by  (5.1)  and  (5.2).   Therefore 
only  internal  gravity  waves  are  present  to  determine  C.   The 
phase  speed  of  the  internal  gravity  waves  is 


C  =  /cpH"  (5.5) 
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where 

g'     is  modified  acceleration  of  gravity  and 
given  by  the  relation 

g'  =  g  A£  (5.6) 

yo 

where 

g     acceleration  of  gravity 

Ap    a  typical  value  for  the  vertical  density 
differences 

p     reference  density 
o 

—2  —3      —3 

Assuming  g  =  1000  cm  sec   ,  Ap  =  10  x  10   gm  cm   and 

—3  ~  —2 

p   =  1.012  gm  cm   gives  a  value  for  g'  =  9.88  cm  sec 

Integrations  are  carried  out  with  an  8— minute  time  step  for 

a  grid  spacing  of  5  km,   without  any  stability  problem. 

B.   THE  FINITE  DIFFERENCE  EQUATIONS 

The  numerical  method  is  set  down  in  terms  of  the  non- 
dimensional  variables.   The  finite  difference  scheme  is  based 
on  a  three— dimensional  array  of  points  with  indices  i,j,k. 
Time  step  is  indicated  by  (n)  as  a  superscript.   In  the  nota- 
tion the  letters  i,j,k  are  always  integers. 

Horizontal  spacing  is  such  that  Ax  =  10  km  and  Ay  =  5  km, 
and  uniform  in  the  x  and  y  directions.   The  horizontal  grid 
pattern  is  shown  in  figure  9.  Where  horizontal  velocity  com- 
ponents (u,v)  are  defined  at  circled  points,  vertical  veloc- 
ity (w)  and  sigma— t  (a)  are  defined  at  dot  points. 
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J:JM 


H  A.  fr-j, 


J:1 
1:1 


T 


o  «.»  Point 
Point 


l:IM 


Figure  9.   Placement  of  variables  on 
horizontal  grid  plane.   The  open  cir- 
cles denote  the  definition  points  of 
the  horizontal  velocity  components 
(u,v)  and  dots  denote  definition  points 
of  (a,w) .   The  distance  between  adja- 
cent grid  point  is  Ax  =  10  km  and 
Ay  =  5  km  in  x  and  y  direction,  res- 
pectively.  IM  and  JM  have  values  21 
and  11  respectively. 
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For  an  arbitrary  interior  point,  where  (u,v)  are  stored 
longitude  (x)  and  latitude  (y)  are  defined  by 


x±   +  \   =  ~  [1  +  (i-1)]       i  =  1,2,...,  IM-1 


yj  +  \   =  ^  [1  +  (j-1)]       j  =  1,2,...,  JM-1 


(5.7) 


The  vertical  pattern  of  variables  is  shown  in  figure  10. 
The  vertical  index  indicated  with  k  and  the  vertical  velocity 
w  are  given  along  the  surface  and  bottom  and  along  the  layer 
boundaries.   The  variables  (u,v,a)  are  located  at  the  kth 
level.   Depth  of  each  level  is  defined  as 
Az, 

zi  -  "  TT 

(5.8) 


m=k 

Z,  = 

j. 

m=2 


k  =  Z1  -   £   <Azm  +  Azm-1)/2   k  =  2/3... KM 


Layer  thicknesses  for  each  level  are  given  in  table  3. 
Table  3.   LAYER  THICKNESSES 


AZ, 

-^    0.25   0.25  0.3125  0.6875  1.0    1.5    6.0   15.0    25.0 


To  resolve  the  surface  layer  more  accurately,  four  levels 
are  located  within  2  0  meters.   Close  to  the  bottom  a  coarser 
grid  is  used.   Resolution  in  the  surface  layer  is  greater  than 
in  the  bottom  layer. 
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1^_ _^__ 0 

1        : )L±1 2.5 

w 

2        !Jli! 7.5 

w 

3      \hl* 12.5 

w 

4     ".v.y 20.0 

w 

5     iLIiI 40. 0 

w 

6      »£! 6  0.0 

6^  —     ___.!*___     __ 

K-7      «!*£ I0  00 


71/2_    _     _    _     «- 


8      !i^ll 3  0  00 


w 


KMS9      'L^l 700.0 


9\ 


* .     10  0  0.0 

/iriniim,w/iiin)ii,'m,wniiii/i/i)iii!iii)ii/iii/iii)ih!::),iin,n<       w 


Figure  10.   Vertical  structure  of  the 
model  (u,v,a)  are  defined  at  integer  K 
values  and  w  is  defined  along  the  layer 
boundaries . 
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The  finite  difference  form  of  the  hydrostatic  equation 
has  the  form 

pijl  "  (0"ao>i,j,35  A\ 

Pijk  =   Pij,k-1  +  f^o'ijk-Js  -L\-k    k=2'3'KM        <5-9» 

and  the  vertical  average  of  the  departure  pressure- is 
calculated  with  the  relation 


k=KM 

3  .jk  ""k 


P.  .   =  T,     P.  .,  AZ,  (5.10) 


k=l 
The  departure  P ■  is  calculated  from 


P!  ..  =  P.  ..  -  P.  .  (5.11) 

13k    13k    13 


which  is  the  finite  difference  analog  of  (4.23). 

The  finite,  difference  approximation  of  the  continuity 
equation,  in  which  the  vertical  velocity  between  the  levels 
is  calculated,  has  the  form 

Wijk+^  -  "ijk-!,  +  [UX  +  Vijk  AZk  (5-12) 

where 


U      .  .,     =  ^(u      .    ,     +  u      ._,).■, 
x   13k  x   j+%  x   3— *g   lk 


V      .  ..    =  3g(v     i+h  +  vy  i-Vik 

v   ilk  v  *  J 


(5.13) 


y  ID 
and 


ux   ij+3gk  =:   Ax{ui+h        u±-h)j+hk 


vy   i+?sjk        hy{vj+k        Vj-Vi+Jsk 
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(5.14) 


which  are  analogs  of  (-^-)  and  (•?— )  respectively.   Gradients 

dx       dy 

of  horizontal  velocity  components  u  and  v   are  defined  at  half- 
integer  grid  points  west  and  south  of  the  u  and  v  storage 
points  shown  in  figure  11a. 

Integration  from  the  surface  downward,  using  surface 
boundary  condition  w=0,  is  done  by 


r.  ..  _,,  =  y  (U  +  V  )  . . 

i jk+*s    La         x    y  lj 


(5.15) 


The  formula  for  the  finite  difference  approximation  of 
the  equation  (4.19)  by  which  u  is  predicted  is 


n+1    n-1  n    P'    .._  -  P!  ... 

ru     —  u    i  1   r  1+1  j+l 1  j+1 

1     2At     Ji+3gj+%k  '   '  2Ro  L        Ax 

+  Pi+1  J  ""  PiJin  .  ^H  rux  i+1  j+^  ""  Ux  i  j+^ 
Ax     Jk   Ro  L        Ay 

Uy  i+^s  j  +  1  ~  Uy  i+^j.n-1   ^v  ruz  k-1^  ""  uz  k+^.n-1 


Ay  Jk     Ro  L     AZ.  i+^  j  + 


"2 


+  55  [cl  v"+1  +  C2  v"  ^i+a,  j+%k  (5-16) 

whe"    "yj+l=  J   'Ay  **>!+*  <5'17' 

Equations  (5.16)  and  (5.17)  represent  (-A  and  (yi)  and 
are  located  half  grid  distance  south  and  above  the  velocity 
(u)  storage  points,  respectively,  as  shown  in  figure  11  a— b. 
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The  formula  for  the  finite  difference  approximation  of 
the  equation  (4.20)  in  which  v  is  predicted  has  the  form 


rVn+l  _  vp-li  j_   .p'i+l   1+1  -  P'i+lj 

1  2At  * ±+hj+kk  2Ro    L  Ay 

v   ij+1  i^n       ~H    ,   x   l+l   j+k  x   13+^ 

Ay  Jk        Ro    l  Ax 


Vy   i+^j+1       Vy  i+^jn  (n-1) 
+  Ay  Jk 


fv    rVz   k-%  "  V: 

Ro    l  AZ,  Ji+^j+^ 


JV    rvz   k-jg  z   k*1^ 


Vr  Lr     y  1     r^,i    n+1    ,    „^   n— 1, 

+  55  h-  t*        -  *5  [c1u        +  c2u       ]i+%j+^k  (5.19) 

r        i+-$ 


where 


v     . . ,  =  Ii±i^ — litk)  J 


x   i+1  Ax 


vk-i  -  vk. 


'j+k 


VZ    kH*    "  AZk_^        )i+3gj+3g  (5'21) 


9v 
Equations  (5.20)  and  (5.21)  are  analogues  of  (-r— )  and 

(^__)  and  are  stored  south  and  above  the  velocity  storage 

dZ 

points,  as  shown  in  figure  11a— b. 
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The  finite  difference  approximation  of  the  sigma— t  equa- 
tion has  the  following  form, 

n+ 1    n— 1  ,    .        -.     o      ■  .  i  ■  .  i  ™"  o    1,1 

r£_-Z_JL_i    =  -a    (n)+  1  /  x  x+ki+k       x  i-kj+k 
1  JEt  Jijk  '    ljk     Pe^         Ax 


°y  i+ki+k     °Y   i+^j-^n-1 
Ay         'k 


JL  _L  f    ,qz   k+k       °z   k-k.  n-1 
Kv  Ro     v(  AZk  ; i+kj+k 


+   6c(an+1)ijk  (5-22> 


where 


Aijk  ■  4S[(ui+%j+%  +  ui+M-«i)  •  (aij  +  "i+lj'k 


-   iu±-hi+h  +  ui-^j-^'  •  <°i-lj   +  ^j'k1 


+  4Syt(vi+y+%  +  vi-%j-%>  •  (aij-l  +  aij'k 


-  <vi+!,j-% +  viJ,j-J,>  •  <°ij-l +  "ij'k1 


+  3SzkI(wijW(crijk  +  0ijk-l> 


-  (wijk+^» • <°ijk  +  aijkfl>1  <5-23> 
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and 


ax  i+^jk  "  Ax  (ai+l    ai)jk 


ay  ij+^k  "  Ay"  (aj+l   aj}ik  (5.24) 


(<y~,,  -  cj_^,). 


z  ijk+%   Azk+1    zk    zk+1  ij 


As  shown  in  figure  (11  a— b)  a      .  ,  and  a   ._,  and  a_,_1 
are  defined  at  half  grid  distances  south,  east  and  above  the 
sigma— t  storage  points. 

.  Before  proceeding  to  a  new  time  step  convective  adjustment 
is  applied.   It  is  difficult  to  state  the  convective  adjust- 
ment mechanism.   With  this  mechanism  the  density  structure  is 
forced  to  remain  stable.   It  may  be  expressed  by 


0        aZ  <  ° 

<5   =  O  (5.25) 

c    °°  ~ 

az  >  0 


In  case  of  unstable  stratification,  vertical  mixing  be- 
comes effectively  infinite.   At  each  time  step  this  infinite 
mixing  is  included  in  numerical  solution  by  testing  sigma— t 
profile  for  unstable  lapse  rate  (<?i-  i  >  cr,  )  .   If  this  condi- 
tion exists  new  values  of  sigma— t  for  those  two  layers  are 
set  equal  to  the  vertical  average  sigma— t  of  the  two  layers 
instantaneously.   This  process  is  repeated  until  complete  sta- 
bility is  reached  for  the  entire  layer. 
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C.   SPECIFICATION  OF  BOUNDATION  CONDITIONS 

As  shown  in  Figure  3—1  cr,w  are  variables  which  are  de- 
fined at  lateral  boundaries;  also  horizontal  gradients  of 
(u,v)  are  defined  at  these  lateral  boundaries.   The  boundary 
conditions  are  satisfied  by  the  following  relations  which 
are  given  for  eastern  and  western  boundaries.   Similar  con- 
ditions exist  for  north  and  south  boundaries  with  the  excep- 
tion of  open  boundaries. 

Normal  components  are  set  equal  to  zero  by 


_  2u. 
Ux  i=1     Ax'i=% 


2u> 
Ux  i=IM    Ax'i=IM-% 


(5.26) 


In  the  momentum  equation  zero— slip  condition  is  applied 
to  the  valocity  tangent  to  the  boundary.   On  the  other  hand, 
when  computing  the  advection  term  in  the  density  equation 
free  slip  condition  exists  for  the  velocity  tangent  to  the 
boundary.   This  implies  the  following  conditions 


u.  ,   =   u.  nl 
i=l      i=l% 


i=IM     IM— \ 


(5.27) 
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Definition  of  vertical  gradients  of  (u,v)  allows  the 
writing  of  the  vertical  stress  term  in  flux  form  easily. 
Boundary  conditions  at  the  surface  and  bottom  are  defined 
in  the  following  manner: 

UZ   k  =  1   =  °  (5.28) 

Zonal  wind  stress  is  assumed  zero  and  only  a  constant 

meridional  wind  stress  is  defined. 

,    V  L 

VZ  k  -  2  "  -if   TsY  (5-29) 

On  the  other  hand  it  is  assumed  that  horizontal  velocities 
vanish  at  the  bottom.   The  following  conditions  force  veloc- 
ities (u,v)  to  reach  zero  at  the  bottom 

2ui 

UZ    k=KM+?2  '        Az'k=KM 

(5.30) 

v  =      2V| 

Z    k=KM+%  Az'k=KM 

The  boundary  conditions  of  zero  mass  fluxes  across  the 
side  walls  and  bottom  are  satisfied  by  imposing  the  following 
conditions: 

e™>i=%  =  "(<"»)  1=3/2]  jk 

(5.31) 

<au)i=IM+!5   =   -(OU)i=IM-JS]jk 
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Similar  conditions  exist  at  north  and  south  boundaries. 

In  the  diffusion  term  zero  flux  and  insulation  are 
satisfied  by  the  following  conditions: 


a      =  —  a 
x±=%  xi=3/2 


(5.32) 


a       =  — a 
xi=IM+%    xi=IM-^ 

where  a   is  defined  at  half  grid  distances  west  of  the 
x  - 

sigma— t  storage  points,  as  shown  in  figure  11a.    Similar 
conditions  exist  at  north  and  south  for  meridional  gradients 
of  sigma— t,  with  the  exception  of  open  boundaries  where 
sigma— t  is  defined  outside  the  domain. 

At  the  sea  surface  downward  flux  is  absent. 

°Z   k  =  1   =  °  (5.33) 

The  sigma— t  equation  is  solved  for  only  the  upper  six 
levels.   It  is  assumed  that  sigma— t  for  the  rest  of  the  depth 
range  has  a  value  a    ,  which  is  held  constant  during  integra— 
tion.   This  enters  the  system  in  the  calculation  of  the 
vertical  gradient  of  sigma— t  one  half  grid  distance  below 
the  sixth  level, 


"Z   k=6h     =  ^>k=6  <^4> 


and  o        has  a  value  28.5  in  sigma— t  units. 
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VI.   RESULTS 

Computations  are  carried  out  for  different  values  of  the 
parameters  with  an  integration  period  of  24  hours.   Parameters 
of  these  runs  are  given  in  Table  2.   Results  of  the  model  are 
given  for  the  fields  of  horizontal  velocity  components  (u,v) , 
vertical  velocity  (w)  and  density  (sigma— t)  in  dimensional 
form. 

The  horizontal  velocity  field  at  z  =  — 2.5  m  is  given  in 
Figure  12.   The  horizontal  stress  exerted  on  the  sea  surface 
by  southerly  wind  causes  the  surface  water  to  move  in  a 
generally  North— east  direction.   This  is  a  consequence  of 
movement  of  surface  water  as  an  Ekman  layer,  drifting  to  the 
right  of  the  wind  stress  in  the  Northern  Hemisphere.   Since 
the  curl  of  the  wind  is  zero  divergence  of  the  Ekman  drift  is 
also  zero.   But  horizontal  velocity  is  strongly  convergent  and 
divergent  near  the  boundaries.   Northeastward  drift  takes 
place  over  the  entire  basin  at  this  depth.   The  strongest  flow 
has  a  magnitude  20  cm/sec  in  the  central  part  of  the  basin  and 
diminishes  toward  the  boundaries.   Near  open  boundaries 
velocities  are  smaller  compared  to  the  velocities  at  neighbor- 
ing grid  points.   The  southward  flow  at  these  points  repre- 
sents the  exchange  with  the  Black  Sea  water  at  the  north  and 
Mediterranean  water  at  the  south.   The  magnitudes  of  these 
currents  are  5  cm/sec  and  2  cm/sec  at  grid  points  just  near 
the  entrances  of  the  Straits  of  Canakkale  and  the  Strait  of 
Istanbul.   The  applied  meridional  wind  stress  has  a  value  of 
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—2 

0.3  dyn  cm   which  is  not  strong  enough  to  move  surface  water 

to  the  north  as  it  does  at  neighboring  grid  points. 

In  the  second  layer  z  =  —7.5  m  the  horizontal  velocity 
changes  direction  to  the  right  and  becomes  weaker.   Horizontal 
velocities  at  that  level  are  shown  in  figure  13w;  they  have  a 
magnitude  5  cm  sec   .   Flow  adjacent  to  the  straits  is  still 
southerly  and  has  a  magnitude  5  cm  sec   and  3  cm  sec   at 
the  vicinities  of  the  Straits  of  Istanbul  and  Canakkale 
respectively. 

The  circulation  pattern  at  lower  levels  differs  from  the 
upper  level  flows  both  in  magnitude  and  direction.   These 
differences  can  be  observed  in  the  circulation  pattern  at 
the  z  =  —40.0  m  level  (Figure  14).   The  pattern  is  not  irregu- 
lar compared  to  the  upper  layers. 

The  northeastward  surface  drift  in  the  first  layers  is 
accompanied  by  a  westward  and  partly  southwestward  flow  in 
the  bottom  layers.   A  northward  flow  adjacent  to  the  Strait 
of  Canakkale  brings  Mediterranean  water  into  the  basin.   A 
similar  northward  flow  of  Mediterranean  water  does  not  appear 
at  the  northern  strait  due  to  the  shallow  sill  depth  at  that 
point. 

Flow  in  the  layers  close  to  the  bottom  becomes  weaker. 
The  influence  of  the  circulation  at  the  straits  on  the  general 
pattern  is  small.   The  circulation  pattern  at  a  depth  700.0  m 
is  shown  in  figure  15.   The  maximum  magnitude  is  0.033  cm  sec 

One  advantage  of  the  numerical  modeling  study  is  that 
it  gives  estimates  of  important  oceanographic  variables  that 


60 


\ 

-t 

L 

/ 

cm 

1 

/ 

/ 

CM 

CN 

/ 

c-* 

\ 

— t 

\ 

-1 

I 

1 

I 

1 

PI 

i 

/ 

<M 

/ 

iM 

/ 

\ 

X 

\ 

1 

1 

1 

/ 

/ 

/ 

/ 

\ 

—1 

l-i 

*S 

**> 

•V 

*j 

tO 

CM 

^M 

O 

-t 

\ 

/ 

/ 

/ 

! 

/ 

/ 

L 

L 

\ 

I 

/ 

/ 

/ 

/ 

/ 

/ 

I 

\ 

•^ 

■* 

•» 

Ui 

-«»• 

•<»■ 

en 

<M 

—i 

in 

\jy 

y 

/ 

/ 

/ 

/ 

/ 

/ 

1 

Ui 

ui 

Ui 

ui 

■V 

■♦ 

CM 

—i 

rr- 

tr~ 

/' 

Ui 

/ 

/ 

•<• 

/ 

/ 

I 

—i 

r 

/ 

/ 

cm 

I 

/ 

/ 

/ 

/ 

1 

I 

/ 

1 

/ 

ui 

/ 

Ui 

/ 

/ 

1 

\ 

\ 

—1 

1 

/ 

1 

Ui 

/ 

Ui 

/ 

Ui 

/ 

US 

/ 

L 

\ 

\ 

/ 

i 

1 

/ 

Ui 

/ 

Ui 

/ 

/ 

I 

tT3 

\ 

\ 

/ 

t*5 

i 

/ 

/ 

/ 

/ 

/ 

/ 

\ 

\ 

/ 

i 

7 

/ 

1 

7 

/ 

/ 

\ 

m 

■«• 

ui 

ui 

L<1 

Ui 

•><< 

en 

— i 

V 

—1 

/ 

/ 

/ 

Ui 

/ 

Ui 

/ 

Ui 

/ 

Ui 

/ 

I 

\ 

-^ 

/ 
*1 

/ 

/ 

US 

/ 

/ 

/ 

/ 

/ 

\ 

\ 

-1 

/ 

/ 

/ 

1 

Ui 

t 

t 

1 

i 

I 

\ 

— i 

L 

L 

i 

/ 

/ 

1 

I 

V; 

L 

\ 

\ 

/ 

i 

/ 

! 

/ 

\ 

I 

\ 

^, 

— i 

— ( 

cm 

*i 

pi 

«0 

&i 

CM 

CM 

/• 

\ 

— i 

/ 

i 

/ 

CM 

i 

L 

to 

/ 

—I 

\ 

\ 

cm 

/ 

/ 

CM 

s 

rt 

/ 

CM 

I 

/ 

/ 

*s 

> 


C 
O 
O 
<D 
U) 

U 
O 

CO 

u 
o 

-M 
O 

> 

>i 
+J 

•H 
O 
0 

<-i 
CD 
> 


n3  U 

+J  <1) 

C  +J 

0  CD 

N  6 
•H 

U  ^ 

0  • 

5B  r- 


cn 


u 

•H 

fa 


61 


V 

\ 

\ 

X 

\ 

X 

X. 

T 

X 

*V 

°v 

\ 

\ 

\ 

X 

X 

X 

^ 

\ 

*»> 

c^- 

\ 

\ 

X 

T 

T 

X 

T 

X 

*v 

«>s» 

\ 

X 

r 

r 

r 

r 

r 

X 

\ 

■v 

\ 

r 

r 

r 

r 

r 

r 

T 

\ 

r 

r 

r 

r 

r 

r 

r 

r 

T 

\ 

— * 

r 

r 

r 

r 

r 

r 

r 

r 

\ 

o 

r 

r 

r 

r 

r 

r 

r 

T 

\ 

«*v 

^ 

r 

r 

r 

r 

r 

r 

T 

\ 

\. 

x 

r 

r 

r 

r 

r 

r 

T 

\ 

°v 

\ 

r 

r 

r 

r 

r 

r 

t 

\ 

■\ 

X 

r 

r 

r 

r 

r 

r 

X 

\ 

"\ 

\ 

r 

r 

r 

r 

r 

r 

X 

\ 

\ 

X 

r 

r 

r 

r 

r 

r 

X 

\ 

\ 

X 

r 

r 

r 

r 

r 

T 

X 

X 

\ 

X 

r 

r 

r 

r 

r 

T 

X 

<v 

\. 

X 

T 

T 

r 

T 

T 

T 

\ 

^ 

=v. 

\ 

^ 

\ 

X 

T 

T 

T 

X 

o- 

=- 

°v 

\ 

\ 

X 

T 

^= 

o 

o 

o 

"V. 

\ 

T 

\ 

X 

r 

— = 

o 

— 1 

2-t 

0) 

> 

0) 


Xi 

+J 

•H 
»W 

O 
4-1 

CO 

U 
O 
-P 

o 

> 

>1 
•p 

•H 
O 
O 
rH 

>  — 

M 

iH  <D 
(TJ  -P 
-P    0) 

C    g 

o 

N  O 
•H  • 
H  O 
O  <* 

X  -* 


d) 

U 

•H 
fa 


62 


\ 

i 

1 

\ 

\ 

X 

\ 

1 

X 

-V 

-^ 

N*> 

\, 

X 

"V, 

X 

X 

"V 

->• 

-* 

-— 

~^ 

"No 

\. 

\> 

X 

X 

\, 

~Ss 

— 

— 

X, 

V 

\ 

\ 

X 

X 

\ 

\, 

--» 

-^ 

\» 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

^o 

~^ 

\, 

\ 

\ 

1 

\ 

\ 

X 

\ 

^o 

"">o 

\ 

\ 

\ 

i 

\ 

\ 

V 

\ 

"So 

~^> 

\ 

\ 

1 

1 

\ 

\ 

X 

\ 

Nj 

N, 

\ 

\ 

X 

\ 

\ 

\ 

V. 

\ 

N. 

N> 

\ 

V. 

\ 

\ 

\ 

\ 

\ 

\ 

"--« 

"^ 

\ 

V 

\ 

\ 

\ 

\ 

\ 

V 

N. 

"No 

\ 

\ 

V 

V 

\ 

X 

N. 

--«> 

-« 

~^ 

\ 

V 

V. 

V. 

\ 

\ 

N. 

-° 

y 

^» 

\ 

\ 

X 

\ 

i 

1 

J 

<^- 

o 

"V 

\ 

\ 

\ 

X 

A 

1 

4 

J 

I 

^> 

\ 

\ 

\ 

\ 

\ 

I 

\ 

\ 

~Sq 

■— 

\> 

\ 

\ 

\ 

\ 

\ 

\, 

N» 

-* 

--, 

N» 

V 

X 

\ 

X 

V 

s, 

--X3 

— » 

-O 

N> 

N. 

X 

X 

X 

\> 

\> 

"N> 

■^ 

V 

1 

i 

X 

X 

\ 

\ 

\ 

\ 

N> 

rH 

0) 

> 


+J 
c 

•H 

u 
o 

w 
u 
o 
-p 
o 

> 

>1 
+J 
•H 
O 
0 

rH       • 
Q)  ^ 

>  u 

rH    -P 

t* 

0  O 

N  • 

•H  O 

U  O 

0  I"- 

S  — 


in 


U 

•H 

P4 


63 


cannot  be  obtained  by  measurements.   One  such  variable  is 
the  vertical  velocity.   The  vertical  velocity  at  the  base 
of  the  first  layer  is  given  in  Figure  16.   The  magnitude  of 
the  vertical  velocity  differs  markedly  between  the  boundaries 
and  the  interior.  Rising  motion  (upwelling)  is  taking  place 
in  the  south  and  southwest  and  sinking  motion  (downwelling) 
is  taking  place  in  the  north— northeast.   Vertical  velocity 
is  mainly  determined  by  the  Ekman  drift  current  intersecting 
the  boundaries.   Southerly  wind  stress  produces  upwelling 
and  downwelling  at  the  south  and  north,  respectively.   Up- 
welling extends  to  the  north  on  the  western  side  and  down- 
welling to  the  south  on  the  eastern  side.   The  maximum  magni- 
tude of  the  vertical  motion  is  358  cm  day   .   Strong  upwell- 
ing and  downwelling  take  place  adjacent  to  the  straits.   This 
is  consistent  with  the  nature  of  the  horizontal  current 
system  in  the  vicinity  of  the  straits. 

Vertical  velocity  at  the  bottom  of  the  second  layer, 
z  =  — 7.5  m,  is  shown  in  Figure  17.   The  general  pattern  is 
similar  to  the  vertical  velocity  at  the  base  of  the  first 
layer.   Strong  vertical  velocity  is  present  along  the  lateral 
walls  to  the  southwest  and  northeast  where  it  has  a  magnitude 
+  300  cm  day   .   With  depth  the  effect  of  wind  on  vertical 
velocity  decreases  due  to  the  vertical  stability.   The  verti- 
cal velocity  pattern  at  the  base  of  the  eighth  level  is 
shown  in  Figure  18. 

The  sigma— t  pattern  at  the  first  level,  2.5  m  is  shown 
in  Figure  19.   This  pattern  at  z  =  — 2.5  m  is  specified  by 
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the  surface  boundary  condition.   At  this  level  the  surface 
sigma— t  is  a  function  of  latitude  only  and  held  constant 
during  the  integration  period. 

Deviations  from  the  surface  pattern  begin  to  appear  just 
below  the  first  level  at  7.5m.   The  sigma— t  pattern  for  the 
second  level  is  shown  in  Figure  20.   A  close  relationship 
is  present  between  this  sigma— t  pattern  and  the  vertical 
velocity  field  shown  in  Figure  17.   At  the  north  and  northeast 
sinking  motion  brings  in  low  density  and  water  and  at  the 
south  rising  motion  brings  in  high  density  water  to  this 
depth.   This  pattern  is  more  strongly  tied  to  the  surface 
boundary  conditions  on  sigma— t  than  at  other  levels.   There 
is  a  low  density  water  pool  at  the  north.   As  a  consequence 
of  upwelling  a  high  dense  water  pool  is  present  at  the  south 
extending  in  an  east— west  direction.   Magnitudes  of  sigma— t 
are  16.5  and  18.0  at  north  and  south  respectively. 

The  sigma— t  pattern  at  the  fourth  level  z  =  — 20  m  is 
shown  in  Figure  21.   The  effect  of  vertical  motion  appears 
as  a  low  density  water  pool  at  the  northeast.   Below  this 
level  z  =  — 40   (not  shown) ,   sigma— t  is  uniform  and  has  a 
value  28.5.   This  is  consistent  with  the  sigma— t  pattern 
given  in  Figure  6. 

Further  insight  into  these  results  may  be  gained  by  examin- 
ing the  meridional  and  east— west  cross  sections  of  sigma— t 
and  the  (u,v)  components  of  the  horizontal  velocity.   There 
is  also  the  possibility  of  comparing  the  predicted  sigma— t 
pattern  to  the  observations. 
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Meridional  cross  sections  of  sigma— t  at  different  longi- 
tudes are  shown  in  Figures  22—24.   The  general  patterns  are 
similar.   High  density  water  rises  to  the  surface  and  strati- 
fication becomes  shallower  further  south.   Below  40  m  density 
is  constant.   This  is  consistent  with  the  observed  sigma— t 
field.   At  the  west  and  south  very  strong  upwelling  and  weak 
downwelling  at  north  are  shown  in  Figure  22.   In  the  central 
region  there  is  upwelling  in  the  south  and  downwelling  in  the 
north  as  shown  in  Figure  23.   Further  east  the  effect  of 
strong  downwelling  is  shown  by  the  rate  of  tilting  of  the 
isopycnals  at  the  north  compared  to  the  south  in  Figure  24. 

Meridional  cross  sections  of  zonal  velocity  in  the  east, 
central  and  western  parts  of  the  basin  are  shown  in  Figures 
25—27.   The  zonal  flow  at  the  surface  is  eastward  and  changes 
direction  with  depth.   Zonal  flow  is  westward  at  greater 
depths.   Strong  flow  is  present  at  the  surface  and  decreases 
rapidly  with  depth.   Variation  with  longitude  can  be  observed 
by  comparing  the  figures.   In  the  central  part  zonal  velocity 
has  a  magnitude  16.0  cm/sec  at  the  surface  and  decreases 
both  laterally  and  vertically.   Below  2  0  meters  flow  is  west- 
ward and  has  a  magnitude  0.6  cm  sec 

East— west  cross  sections  of  sigma— t  patterns  for  three 
latitudes  are  shown  in  Figures  28—30.   A  regular  stratifica- 
tion is  deformed  by  upwelling  and  downwelling  of  the  eastern 
and  western  sides  of  the  basin  (Figure  29) .   Departures  are 
present  in  the  north  and  south,  as  shown  in  Figures  28  and  29 
respectively.   Intrusion  of  less  dense  Black  Sea  water  can 
be  identified  in  Figure  30. 
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The  vertical  cross  sections  of  the  meridional  velocity 
which  correspond  to  the  previous  sigma— t  sections,  are  given 
in  Figures  31—33.   There  is  a  northward  flow  at  the  surface 
in  each  section.   The  cross  section  near  the  southern  bound- 
ary shows  a  southerly  flow  below  the  surface  layer.   Flow 
changes  direction  at  about  20  meters.   In  the  vicintiy  of 
the  Strait  of  Istanbul  the  effect  of  the  dpen  boundary  is 
indicated  by  a  weak  northerly  flow.   There  is  a  southerly 
flow  below  that  and  at  depth  about  30  m  another  reversal 
takes  place.   Weak  northward  flow  is  associated  with  the 
southerly  wind.  The  southward  and  northward  flows  at  depth 
are  consistent  with  observations. 
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VII.   CONCLUSIONS 

In  order  to  investigate  the  dynamics  of  the  Sea  of  Mar- 
mara a .circulation  model  is  formulated  based  on  the  hydro- 
static and  Boussinesq  approximations.   Since  very  little  is 
known  about  the  turbulence  characteristics  of  the  ocean  or 
adjacent  seas,  viscosity  and  conductivity  are  replaced  by 
new  terms  representing  the  contributions  of  the  smaller  scale 
motions  to  the  exchange  of  momentum  and  density.   Horizontal 
eddy  transfer  of  momentum  is  accomplished  using  a  constant 
eddy  viscosity  and  diffusivity  coefficient.   For  the  exchange 
of  momentum  and  density  in  the  vertical  a  constant  vertical 
diffusion  coefficient  also  is  used.   The  vertical  mixing  of 
density  is  enhanced  by  the  use  of  an  instantaneous  adjustment 
mechanism  to  eliminate  unstable  lapse  rates.   A  lateral  mix- 
ing coefficient  for  density  is  less  than  that  used  for 
momentum.   These  coefficients  should  be  sufficiently  small 
so  as  not  to  obscure  lateral  transfer  of  the  quantities,  and 
a  non  zero  momentum  coefficient  is  needed  to  satisfy  the 
zero  slip  boundary  condition.   In  this  study  constant  values 
of  the  coefficients  are 

7     2—1 

iL  =  .4  x  10   cm  sec    (density) 


and 


A  7     2—1 

M  =   5  x  10   cm  sec    (momentum] 
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The  vertical  diffusion  coefficient  which  is  difficult  to 

evaluate  either  by  measurement  or  in  principle,  is  assigned 

2   —1 

1.6  and  3.2  cm  sec    for  density  and  momentum  respectively. 

Other  parameters  and  constants  are  given  in  Table  II. 

Further  study  of  thermocline  penetration,  water  exchange, 
vertical  mixing  processes  and  current  gradients  might  provide 
a  better  possibility  for  matching  these  coefficients. 

Although  many  important  factors,  such  as  the  irregular 
bottom  topography,  bottom  friction  and  non— linear  terms  in 
the  equation  of  motion  have  been  omitted,  the  results  from 
the  model  are  generally  encouraging.   More  realistic  results 
may  be  achieved  by  incorporating  these  effects  systematically 
in  further  studies. 

_—  p.   Integration  in  the  model  is   carried  out  with  a  constant 
southerly  wind  stress  which  was  chosen  arbitrarily.   Although 
the  scale  of  the  motion  is  small  compared  to  a  characteristic 
length  scale  for  meteorological  features,  wind  stress  may  be 
significantly  variable  over  the  sea  due  to  the  geographic 
location  of  the  sea.   Also  seasonal  variations  of  the  wind 
may  have  a  considerable  effect  on  the  vertical  mean  part  of 
the  current  which  is  neglected  during  this  study. 

The  model  does  not  separately  predict  temperature  and 
salinity  which  are  sometimes  important  quantities.   Diagnostic 
calculation  of  the  density  with  the  aid  of  predicted  tempera- 
ture and  salinity  fields  can  be  done  with  a  very  simple 
modification  of  the  model.   This  would  require  a  calculation 
of  the  heat  and  salt  fluxes  at  the  surface  and  meridional 
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salt  and  heat  separately  in  the  model.   It  would  be  more 
difficult,  however,  to  define  boundary  conditions  for  these 
quantities  at  the  open  boundaries.   It  is  also  difficult  to 
simulate  seasonal  variations  of  these  quantities  at  the 
open  boundaries.   Once  these  are  defined  properly,  based 
on  accurate  observations,  there  are  no  inherent  difficulties 
in  simulating  the  dynamical  processes  in  the  sea  using  the 
general  principles  involved  in  this  model.   Even  though  this 
model  depends  on  the  particular  hypothesis  used  for  hori— 
zontal  and  vertical  eddy  transport  of  heat,  salt  and  momentum, 
it  is  clear  that  valuable  studies  on  small  scale  water  bodies 
can  be  made. 
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APPENDIX  A 
COMPUTER  PROGRAM  DESCRIPTION 

The  computer  program  is  written  in  FORTRAN  IV  and  was 
used  with  the  IBM  360/67  computer  system  at  the  W.  R.  Church 
Computer  Center,  Naval  Postgraduate  School. 

The  overall  program  is  divided  into  two  basic  subpro- 
grams (1)  the  main  program  and  associated  subroutines  (2)  an 
access  program  which  draws  and  writes  the  results  of  the 
first  program. 

The  main  program  consists  of  nine  subroutines  which 
calculate  different  terms  of  the  equation  of  motion  and  the  . 
density  equation,  pressure,  vertical  velocity  and  changes 
variables  for  next  time  step. 

FORTRAN  IV  symbols  for  the  primary  program  and  a  brief 
description  of  the  subroutines  are  given  below: 

UMI       Zonal  velocity  component  at  n— 1  time  step 
U         Zonal  velocity  component  at  n   time  step 
UA1       Zonal  velocity  component  at  n+1  time  step 
VM1       Meridional  velocity  component  at  n— 1  time  step 
U         Meridional  velocity  component  at  n    time  step 
UA1       Meridional  velocity  component  at  n+1  time  step 
SGMl      Sigma— t  at  n— 1  time  step 
SGMT      Sigma— t  at  n   time  step 
SGMT      Sigma— t  at  n+1  time  step 


89 


ADV  Total  advection  term  in  density  equation 

ADV  Local  rate  of  change  of  sigma— t 

ADV  Pressure 

A  Surface  pressure 

PBAR  Vertical  average  pressure 

W  Vertical  velocity 

UXG  Gradient  of  zonal  velocity  in  x  direction 

UYG  Gradient  of  zonal  velocity  in  y  direction 

UZG  Gradient  of  zonal  velocity  in  vertical 

VXG  Gradient  of  meridional  velocity  in  x  direction 

VYG  Gradient  of  meridional  velocity  in  y  direction 

VZG  Gradient  of  meridional  velocity  in  vertical 

SGX  Gradient  of  sigma— t  in  x  direction 

SGY  Gradient  of  sigma— t  in  y  direction 

SGZ  Gradient  of  sigma— t  in  vertical 

TX  Meridional  wind  stress 

TY  Zonal  wind  stress 

DZ  Layer  depth 

DZ  Layer 

BB  Depth  dependent  velocity  at  the  Northern 
fictitious  boundary 

CA  Depth  dependent  velocity  at  the  Southern 

fictitious  boundary 

fictitious  boundary 

BO  Sigma— t  defined  outside  the  domain  at  North 

CD  Sigma— t  defined  outside  the  domain  at  South 

DX  Horizontal  grid  spacing  in  x  direction 

DY  Horizontal  grid  spacing  in  y  direction 

RO  Rossby  number 
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EH 

EV 

PE1 

Wl 

AK 

AKV 

RL 

HH 

DT 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 

SUBROUTINE 
SUBROUTINE 


Horizontal  Ekman  number 
Vertical  Ekman  number 
P^clet  number 
Relative  velocity 
Vertical  eddy  diffusivity 
Vertical  eddy  viscosity 
Curl  factor 
Depth  of  the  basin 
Time  step 
PRES     Calculates  pressure 

Calculates  vertical  velocity 


VERW 
ADVEC 

HDGRD 
DIFFO 


SUBROUTINE  SIGEQ 


SUBROUTINE  HVGRD 


SUBROUTINE  HORV 


Calculates  advection  term  in  the 
density  equation 

Calculates  gradients  of  the  sigma— t 

Calculates  local  time  rate  of  the 
change  sigma— t  by  subtracting  advec- 
tion term  from  the  diffusion  term 

Calculates  new  sigma— t  and  makes 
convective  adjustment  for  unstable 
lapse  rates. 

Calculates  gradients  of  the  hori- 
zontal velocities 

Calculates  horizontal  velocity 
components,  u,v. 


A  descriptive  flow  diagram  of  the  program  is  shown  in 
Figure  34  . 
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(  start  ) 


Layers 
Depth 


Calculate 

Layer 
Thickness 


Meridional  Wind 
Stress    (TY) 


Define  Zonal  Wind 
Stress    (TY) 


Define  Surface 
Pressure 


Define  Parameters 
of  the  Model 


Figure   34.      Descriptive   flow  diagram  of   the  program, 
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Set  Initial 
Velocities  (U,V,W) 


Define  Initial 
Sigma-T  Distribution 


Define  Sigma-T 
at  the  Fictitious 
Grid  Points 


N  «-l 


NMAX 


-0 


U,V,W,SGMT 
(      STOP 


Figure   34.      Continued 


Yes 


SAGLA 


HDGRD 


DZ,BO,CD,BB,CA,SGMT,SGX,SGY,SGZ 


Yes 


HDGRD 


,CA,SGM1,SGX 
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Figure   34.    Continued 
94 


Yes 


No 


ADV 


CALCULATE  VERTICAL 
AVERAGE  PRESSURE  (PBAR) 


ADV(I,J,K)  =  ADV(I,J,K)  -  PBAR 


Figure  34.   Continued 
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HVGRD 


. ..,TY,    UM1,    VMl,    UXG,.. 


HORV 


D4T,DZ,ADV,    UXG , VXG , UYG , VYG ,    UZG,VZG,TX, TY,    U,V,UA1,VA1 
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No    X 


Yes 


HORV 


D2T, . . . ,UM1,VM1,UA1,VA1 


Figure   34.      Continued 
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0 


HORV 


DT, . . . , UM1,VM1,UA1, VA1 


CVAR 


UA1,  U,  VA1,  V,  SGP1,  SGMT,  SGPl,  UM1,  U,  VMl,  V,  SGM1,  SGMT 

I 

w    

CVAR 


U,  UA1,  V,  VA1,  SGMT,  SGP1,  UM1,  U,  VMl,  V,  SGMl,  SGMT 


Figure  34.   Continued 
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Figure  34.   Continued 
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